## Courting the President replication
##
## non-contender analysis 
## Estimand: ATT 
## Dependent variable: presIdeoVote 
##
## Reproduce: 
##      - Table 1 (bottom panel) 
##      - Figure A1 (TASMD & KS statistic)
##
## Kevin Quinn and Guoer Liu
## 11/19/2024

# setwd("") # set at your current working directory
set.seed(91800198)


source("func.R")
library(haven)
library(WeightIt)
library(Matching)
library(cem)
library(data.table)
library(RColorBrewer)
library(lattice)



mydata <- read_dta("NonContender.dta")
mydata <- as.data.frame(mydata)

## drop casepub b/c no variation in noncontender data
mydata <- mydata[, c("presIdeoVote", "nSLtreatFinal0",
                     "judgeJCS", "presDist",
                     "panelDistJCS", "circmed", "sctmed", "coarevtc",
                     "judge")]
mydata <- na.omit(mydata)

cov.names.sub <- c("judgeJCS", "presDist", "panelDistJCS", "circmed",
                   "sctmed", "coarevtc")

treat.form <- "nSLtreatFinal0 ~ judgeJCS + presDist + panelDistJCS + circmed + sctmed + coarevtc"


reg.general.form <- "presIdeoVote ~ judgeJCS + presDist + panelDistJCS + circmed + sctmed + coarevtc"




estimator <- NULL
estimate <- NULL
SE <- NULL
Eff.n <- NULL
Eff.n.ratio <- NULL
DFBETA.min <- NULL
DFBETA.max <- NULL
DFBETA.min.who <- NULL
DFBETA.max.who <- NULL
extrap.treat <- NULL
extrap.control <- NULL




## #############################################################
## general regression
estimator <- c(estimator, "reg.general")

mydata.1 <- mydata[mydata$nSLtreatFinal0==1,]
mydata.0 <- mydata[mydata$nSLtreatFinal0==0,]

lm.out.0 <- lm(reg.general.form, data=mydata.0)

m0 <- predict(lm.out.0, newdata=mydata.1)

estimate <- c(estimate, mean(mydata.1$presIdeoVote - m0))



## SE via nonparametric bootstrap
B <- 5000
n <- nrow(mydata)
estim.bs <- rep(NA, B)
for (b in 1:B){
    bs.inds <- sample(1:n, n, replace=TRUE)
    mydata.bs <- mydata[bs.inds,]

    mydata.bs.1 <- mydata.bs[mydata.bs$nSLtreatFinal0==1,]
    mydata.bs.0 <- mydata.bs[mydata.bs$nSLtreatFinal0==0,]

    lm.out.0 <- lm(reg.general.form, data=mydata.bs.0)
    
    m0 <- predict(lm.out.0, newdata=mydata.bs.1)

    estim.bs[b] <- mean(mydata.bs.1$presIdeoVote - m0)
}
SE <- c(SE, sd(estim.bs))


## get the WATE weights
X1 <- as.matrix(mydata[mydata$nSLtreatFinal0==1, cov.names.sub])
X1 <- cbind(1, X1)
X0 <- as.matrix(mydata[mydata$nSLtreatFinal0==0, cov.names.sub])
X0 <- cbind(1, X0)
n1 <- nrow(X1)
n0 <- nrow(X0)
n <- n1 + n0
treat.vec <- c(rep(1, n1), rep(0, n0))
y1 <- mydata[mydata$nSLtreatFinal0==1, "presIdeoVote"]
y0 <- mydata[mydata$nSLtreatFinal0==0, "presIdeoVote"]
X <- rbind(X1, X0)
y <- as.matrix(c(y1, y0), nrow=n, ncol=1)
ones.n <- matrix(1, nrow=n1, ncol=1)

w1 <- rep(1, n1)
w0 <- t(ones.n) %*% X1 %*% solve(t(X0) %*% X0) %*% t(X0)
w <- as.vector(c(w1, w0))

n.eff.0 <- (sum(w0)^2) / sum(w0^2)
n.eff.1 <- (sum(w1)^2) / sum(w1^2)
n.eff <- (n.eff.1 * n.eff.0) / (n.eff.1 + n.eff.0)
n.eff.max <- (n1 * n0) / (n1 + n0)
Eff.n <- c(Eff.n, n.eff)
Eff.n.ratio <- c(Eff.n.ratio, n.eff/n.eff.max)

DFBETA.out <- DFBETA.weighted(y, treat.vec, w)
DFBETA.min <- c(DFBETA.min, min(DFBETA.out))
DFBETA.max <- c(DFBETA.max, max(DFBETA.out))

min.ind <- which.min(DFBETA.out)
max.ind <- which.max(DFBETA.out)

mydata.noncontender.reg.att <- rbind(mydata.1, mydata.0)
mydata.noncontender.reg.att$w <- w

DFBETA.min.who <- c(DFBETA.min.who, 
                    mydata.noncontender.reg.att[min.ind, "judge"])
DFBETA.max.who <- c(DFBETA.max.who, 
                    mydata.noncontender.reg.att[max.ind, "judge"])

extrap.1 <- sum(abs(w1) * (w1 < 0)) / sum(w1 * (w1 >= 0))
extrap.0 <- sum(abs(w0) * (w0 < 0)) / sum(w0 * (w0 >= 0))

extrap.treat <- c(extrap.treat, extrap.1)
extrap.control <- c(extrap.control, extrap.0)



## END general regression (two regression models)
## #############################################################









## #############################################################
## general regression (constant effect)
estimator <- c(estimator, "reg.constant")
reg.cons.form <- "presIdeoVote ~ nSLtreatFinal0 + judgeJCS + presDist + panelDistJCS + circmed + sctmed + coarevtc"
lm.out.cons <- lm(reg.cons.form, data=mydata)

estimate <- c(estimate, summary(lm.out.cons)$coefficient[2, 1])

## SE via nonparametric bootstrap
B <- 5000
n <- nrow(mydata)
estim.bs <- rep(NA, B)
for (b in 1:B){
    bs.inds <- sample(1:n, n, replace=TRUE)
    mydata.bs <- mydata[bs.inds,]

    lm.out.cons <- lm(reg.cons.form, data=mydata.bs)

    estim.bs[b] <- summary(lm.out.cons)$coefficient[2, 1]
}
SE <- c(SE, sd(estim.bs))


## get the WATE weights
treatment <- mydata$nSLtreatFinal0
y <- mydata$presIdeoVote
n <- length(treatment)
n1 <- sum(treatment == 1)
n0 <- sum(treatment == 0)

X <- as.matrix(mydata[, cov.names.sub])
X.good.cols <- qr(X)$pivot[seq_len(qr(X)$rank)]
ordvec <- order(treatment)
treat <- treatment[ordvec]
X.cov <- X[ordvec, X.good.cols]
X <- cbind(1, treat, X.cov)
X0 <- cbind(1, 0, X.cov)
X1 <- cbind(1, 1, X.cov)

ones <- matrix(1, nrow = n, ncol = 1)
    
w1.full <- t(ones) %*% X1 %*% solve(crossprod(X)) %*% t(X)
w0.full <- t(ones) %*% X0 %*% solve(crossprod(X)) %*% t(X)
w1 <- w1.full[(n0+1):n] - w0.full[(n0+1):n]
w0 <- w0.full[1:n0] - w1.full[1:n0]
w.ordvec <- c(w0, w1) # stacked w
w <- w.ordvec[order(ordvec)] # w in the original order

w1 <- w[treatment == 1]
w0 <- w[treatment == 0]

n.eff.0 <- (sum(w0)^2) / sum(w0^2)
n.eff.1 <- (sum(w1)^2) / sum(w1^2)
n.eff <- (n.eff.1 * n.eff.0) / (n.eff.1 + n.eff.0)
n.eff.max <- (n1 * n0) / (n1 + n0)
Eff.n <- c(Eff.n, n.eff)
Eff.n.ratio <- c(Eff.n.ratio, n.eff/n.eff.max)

DFBETA.out <- DFBETA.weighted(y, treatment, w)
DFBETA.min <- c(DFBETA.min, min(DFBETA.out))
DFBETA.max <- c(DFBETA.max, max(DFBETA.out))

min.ind <- which.min(DFBETA.out)
max.ind <- which.max(DFBETA.out)

DFBETA.min.who <- c(DFBETA.min.who, mydata[min.ind, "judge"])
DFBETA.max.who <- c(DFBETA.max.who, mydata[max.ind, "judge"])

extrap.1 <- sum(abs(w1) * (w1 < 0)) / sum(w1 * (w1 >= 0))
extrap.0 <- sum(abs(w0) * (w0 < 0)) / sum(w0 * (w0 >= 0))

extrap.treat <- c(extrap.treat, extrap.1)
extrap.control <- c(extrap.control, extrap.0)

mydata.noncontender.reg_cons.att <- cbind(mydata, w)

## END general regression (constant effect)
## #############################################################










## #############################################################
## propensity score weighting (logit)
estimator <- c(estimator, "sipw.logit")

glm.out <- glm(treat.form, data=mydata, family=binomial(logit))
mydata$pscores <- glm.out$fitted
mydata$w <- 1 
mydata$w[mydata$nSLtreatFinal0==0] <-
    mydata$pscores[mydata$nSLtreatFinal0==0] /
    (1 - mydata$pscores[mydata$nSLtreatFinal0==0])

ATT.hat <- sum(mydata$w * mydata$nSLtreatFinal0 * mydata$presIdeoVote) /
    sum(mydata$w * mydata$nSLtreatFinal0) -
    sum(mydata$w * (1-mydata$nSLtreatFinal0) * mydata$presIdeoVote) /
    sum(mydata$w * (1-mydata$nSLtreatFinal0))

estimate <- c(estimate, ATT.hat)


## SE via nonparametric bootstrap
B <- 5000
n <- nrow(mydata)
estim.bs <- rep(NA, B)
for (b in 1:B){
    bs.inds <- sample(1:n, n, replace=TRUE)
    mydata.bs <- mydata[bs.inds,]

    glm.out.bs <- glm(treat.form, data=mydata.bs,
                      family=binomial(logit))
    mydata.bs$pscores <- glm.out.bs$fitted
    mydata.bs$w <- 1 
    mydata.bs$w[mydata.bs$nSLtreatFinal0==0] <-
        mydata.bs$pscores[mydata.bs$nSLtreatFinal0==0] /
        (1 - mydata.bs$pscores[mydata.bs$nSLtreatFinal0==0]) 
     
    ATT.hat.bs <- sum(mydata.bs$w * mydata.bs$nSLtreatFinal0 *
                      mydata.bs$presIdeoVote) /
        sum(mydata.bs$w * mydata.bs$nSLtreatFinal0) -
        sum(mydata.bs$w * (1-mydata.bs$nSLtreatFinal0) *
            mydata.bs$presIdeoVote) /
        sum(mydata.bs$w * (1-mydata.bs$nSLtreatFinal0))
    
    
    estim.bs[b] <- ATT.hat.bs
}
SE <- c(SE, sd(estim.bs))


w0 <- mydata$w[mydata$nSLtreatFinal0==0]
w1 <- mydata$w[mydata$nSLtreatFinal0==1]
n1 <- sum(mydata$nSLtreatFinal0)
n <- nrow(mydata)
n0 <- n - n1

n.eff.0 <- (sum(w0)^2) / sum(w0^2)
n.eff.1 <- (sum(w1)^2) / sum(w1^2)
n.eff <- (n.eff.1 * n.eff.0) / (n.eff.1 + n.eff.0)
n.eff.max <- (n1 * n0) / (n1 + n0)
Eff.n <- c(Eff.n, n.eff)
Eff.n.ratio <- c(Eff.n.ratio, n.eff/n.eff.max)

DFBETA.out <- DFBETA.weighted(mydata$presIdeoVote, mydata$nSLtreatFinal0,
                              mydata$w)
DFBETA.min <- c(DFBETA.min, min(DFBETA.out))
DFBETA.max <- c(DFBETA.max, max(DFBETA.out))

min.ind <- which.min(DFBETA.out)
max.ind <- which.max(DFBETA.out)
DFBETA.min.who <- c(DFBETA.min.who, mydata[min.ind, "judge"])

DFBETA.max.who <- c(DFBETA.max.who, mydata[max.ind, "judge"])


extrap.treat <- c(extrap.treat, NA)
extrap.control <- c(extrap.control, NA)


mydata.noncontender.sipw.logit.att <- mydata


## END propensity score weighting (logit)
## #############################################################







## #############################################################
## entropy balancing

estimator <- c(estimator, "entropy.balancing")

eb.out <- weightit(as.formula(treat.form), data=mydata,
                   method="ebal",
                   estimand="ATT")


mydata$w <- eb.out$weights

wls.out <- lm(presIdeoVote ~ nSLtreatFinal0, weights=w, data=mydata)

ATT.hat <- coef(wls.out)[2]

estimate <- c(estimate, ATT.hat)


SE <- c(SE, sqrt(vcov(wls.out)[2,2]))

w0 <- mydata$w[mydata$nSLtreatFinal0==0]
w1 <- mydata$w[mydata$nSLtreatFinal0==1]
n1 <- sum(mydata$nSLtreatFinal0)
n <- nrow(mydata)
n0 <- n - n1

n.eff.0 <- (sum(w0)^2) / sum(w0^2)
n.eff.1 <- (sum(w1)^2) / sum(w1^2)
n.eff <- (n.eff.1 * n.eff.0) / (n.eff.1 + n.eff.0)
n.eff.max <- (n1 * n0) / (n1 + n0)
Eff.n <- c(Eff.n, n.eff)
Eff.n.ratio <- c(Eff.n.ratio, n.eff/n.eff.max)

DFBETA.out <- DFBETA.weighted(mydata$presIdeoVote, mydata$nSLtreatFinal0,
                              mydata$w)
DFBETA.min <- c(DFBETA.min, min(DFBETA.out))
DFBETA.max <- c(DFBETA.max, max(DFBETA.out))

min.ind <- which.min(DFBETA.out)
max.ind <- which.max(DFBETA.out)
DFBETA.min.who <- c(DFBETA.min.who, mydata[min.ind, "judge"])

DFBETA.max.who <- c(DFBETA.max.who, mydata[max.ind, "judge"])


extrap.treat <- c(extrap.treat, NA)
extrap.control <- c(extrap.control, NA)

mydata.noncontender.entropy.balancing.att <- mydata


## END entropy balancing
## #############################################################








## #############################################################
## 1:1 NN matching
estimator <- c(estimator, "NN.matching")

treat.vec <- mydata$nSLtreatFinal0
X <- as.matrix(mydata[, c(cov.names.sub)])
y <- mydata$presIdeoVote

# Matching by linearized propensity score
glm.out <- glm(treat.vec ~ X, family=binomial(logit))
l.pscores <- predict(glm.out)

match.out <- Match(Y = y, Tr=treat.vec, X=l.pscores, estimand="ATT",
                   M=1, ties=FALSE, replace=TRUE, BiasAdjust=FALSE)

estimate <- c(estimate, match.out$est)
SE <- c(SE, match.out$se.standard)

mydata$w <- 0
for (i in 1:nrow(mydata)){
     mydata$w[i] <- sum(match.out$index.control == i) +
            sum(match.out$index.treated == i)
}

w0 <- mydata$w[mydata$nSLtreatFinal0==0]
w1 <- mydata$w[mydata$nSLtreatFinal0==1]
n1 <- sum(mydata$nSLtreatFinal0)
n <- nrow(mydata)
n0 <- n - n1

n.eff.0 <- (sum(w0)^2) / sum(w0^2)
n.eff.1 <- (sum(w1)^2) / sum(w1^2)
n.eff <- (n.eff.1 * n.eff.0) / (n.eff.1 + n.eff.0)
n.eff.max <- (n1 * n0) / (n1 + n0)
Eff.n <- c(Eff.n, n.eff)
Eff.n.ratio <- c(Eff.n.ratio, n.eff/n.eff.max)

DFBETA.out <- DFBETA.weighted(mydata$presIdeoVote, mydata$nSLtreatFinal0,
                              mydata$w)
DFBETA.min <- c(DFBETA.min, min(DFBETA.out))
DFBETA.max <- c(DFBETA.max, max(DFBETA.out))

min.ind <- which.min(DFBETA.out)
max.ind <- which.max(DFBETA.out)
DFBETA.min.who <- c(DFBETA.min.who, mydata[min.ind, "judge"])

DFBETA.max.who <- c(DFBETA.max.who, mydata[max.ind, "judge"])


extrap.treat <- c(extrap.treat, NA)
extrap.control <- c(extrap.control, NA)

mydata.noncontender.NN.matching.att <- mydata

## END 1:1 NN matching
## #############################################################













## #############################################################
## coarsened exact matching
estimator <- c(estimator, "CEM")

cem.out <- cem(treatment="nSLtreatFinal0", data=mydata,
               drop=c("presIdeoVote", "nSLtreatFinal0", "judge",
                      "pscores", "w"))
att.out <- att(cem.out, formula = presIdeoVote ~ nSLtreatFinal0,
               data=mydata)
s <- summary(att.out)


estimate <- c(estimate, s[2,1])
SE <- c(SE, s[2,2])

mydata$w <- cem.out$w

w0 <- mydata$w[mydata$nSLtreatFinal0==0]
w1 <- mydata$w[mydata$nSLtreatFinal0==1]
n1 <- sum(mydata$nSLtreatFinal0)
n <- nrow(mydata)
n0 <- n - n1

n.eff.0 <- (sum(w0)^2) / sum(w0^2)
n.eff.1 <- (sum(w1)^2) / sum(w1^2)
n.eff <- (n.eff.1 * n.eff.0) / (n.eff.1 + n.eff.0)
n.eff.max <- (n1 * n0) / (n1 + n0)
Eff.n <- c(Eff.n, n.eff)
Eff.n.ratio <- c(Eff.n.ratio, n.eff/n.eff.max)

DFBETA.out <- DFBETA.weighted(mydata$presIdeoVote, mydata$nSLtreatFinal0,
                              mydata$w)
DFBETA.min <- c(DFBETA.min, min(DFBETA.out))
DFBETA.max <- c(DFBETA.max, max(DFBETA.out))

min.ind <- which.min(DFBETA.out)
max.ind <- which.max(DFBETA.out)
DFBETA.min.who <- c(DFBETA.min.who, mydata[min.ind, "judge"])

DFBETA.max.who <- c(DFBETA.max.who, mydata[max.ind, "judge"])


extrap.treat <- c(extrap.treat, NA)
extrap.control <- c(extrap.control, NA)

mydata.noncontender.CEM.att <- mydata

## END coarsened exact matching
## #############################################################


## #############################################################
## TABLE

## Table 1 (bottom panel) ###


## noncontender

noncontender.att <- data.frame(estimator=estimator,
                           estimate=estimate,
                           SE=SE,
                           #Z = estimate / SE,
                           Eff.n=Eff.n,
                           Eff.n.ratio=Eff.n.ratio,
                           DFBETA.min=DFBETA.min,
                           DFBETA.min.who=DFBETA.min.who,
                           DFBETA.max=DFBETA.max,
                           DFBETA.max.who=DFBETA.max.who,
                           extrap.control=extrap.control,
                           extrap.treat=extrap.treat)
print(noncontender.att)
#library(xtable)
#xtable(noncontender.att, digits=c(0, 0, rep(3, 5), 0, 3, 0, rep(3, 2)))


## #############################################################
## FIGURES 

## Figure A1: TASMD (left panel) ###

mydata.raw <- mydata
mydata.raw$w <- 1

data.list <- list(mydata.raw,
                  mydata.noncontender.reg.att,
                  mydata.noncontender.reg_cons.att,
                  mydata.noncontender.sipw.logit.att,
                  mydata.noncontender.entropy.balancing.att,
                  mydata.noncontender.NN.matching.att,
                  mydata.noncontender.CEM.att)

TASMD.noncontender.data <- as.data.frame(matrix(nrow = length(cov.names.sub), 
                                                ncol = length(data.list)))
colnames.var <- c("Raw Data",
                  "Regression (MRI)",
                  "Regression (URI)",
                  "PScores",
                  "Ebal",
                  "NNmatch", 
                  "CEM")

colnames(TASMD.noncontender.data) <- colnames.var
rownames(TASMD.noncontender.data) <- cov.names.sub

for(i in 1:length(data.list)){
    mydata <- data.list[[i]]
    for (xvar in cov.names.sub){
        tasmd <- TASMD(x=mydata[,xvar], w=mydata$w,
                       treatment=mydata$nSLtreatFinal0)
        TASMD.noncontender.data[xvar, i] <- tasmd
    }
}

xvar.names <- c("JCS score", "Ideo. Distance", "Panel JCS",
                 "Circuit med", "SC med", "Court reversal")
TASMD.noncontender.data$xvar <- xvar.names


TASMD.noncontender.long <- melt(setDT(TASMD.noncontender.data),
                            id.vars="xvar", variable.name = "method")


colors <- brewer.pal(7, "Dark2")
colors <- replace(colors, c(1, 2), colors[c(2, 1)])

#pdf("TASMD-noncontenders.pdf", height=7, width=6, family = "Times")
trell.par <- trellis.par.get() ## get default settings
trellis.par.set(superpose.symbol =
                    list(pch = c(16, 0:3, 18, 4), cex = 1,
                         col = colors                                    
                         ))
print(dotplot(xvar ~ value, 
              data = TASMD.noncontender.long,
              groups = method,
              xlab="TASMD",
              xlim = c(-0.1, 0.9), #@ same scale as the contender TASMD
              scales=list(y=list(cex=1.2), x=list(cex=1.2)),
              panel=function(x, y, ...){
                  panel.abline(h=1:42,col="grey95", ...)
                  panel.abline(v=0.1, col="grey85")
                  panel.xyplot(x, y, ...)
              },
              key=simpleKey(levels(TASMD.noncontender.long$method),
                            space="right")
              ))
trellis.par.set(trell.par)  ## go back to default
#dev.off()


### Figure A1: KS statistics (right panel) ###

KS.noncontender.data <- as.data.frame(matrix(nrow = length(cov.names.sub), 
                                                ncol = length(data.list)))

colnames(KS.noncontender.data) <- colnames.var
rownames(KS.noncontender.data) <- cov.names.sub

for(i in 1:length(data.list)){
    mydata <- data.list[[i]]
    for (xvar in cov.names.sub){
        KS <- wtd.ks(x=mydata[,xvar], w=mydata$w,
                       treatment=mydata$nSLtreatFinal0)
        KS.noncontender.data[xvar, i] <- KS
    }
}

KS.noncontender.data$xvar <- xvar.names


KS.noncontender.long <- melt(setDT(KS.noncontender.data),
                            id.vars="xvar", variable.name = "method")


colors <- brewer.pal(7, "Dark2")
colors <- replace(colors, c(1, 2), colors[c(2, 1)])

#pdf("KS-noncontenders.pdf", height=7, width=6, family = "Times")
trell.par <- trellis.par.get() ## get default settings
trellis.par.set(superpose.symbol =
                    list(pch = c(16, 0:3, 18, 4), cex = 1,
                         col = colors                                    
                         ))
print(dotplot(xvar ~ value, 
              data = KS.noncontender.long,
              groups = method,
              xlab="KS Statistics",
              xlim = c(-0.05, 0.49), #@ same scale as the contender KS
              scales=list(y=list(cex=1.2), x=list(cex=1.2)),
              panel=function(x, y, ...){
                  panel.abline(h=1:42,col="grey95", ...)
                  panel.abline(v=0, col="grey85")
                  panel.xyplot(x, y, ...)
              },
              key=simpleKey(levels(KS.noncontender.long$method),
                            space="right")
              ))
trellis.par.set(trell.par)  ## go back to default
#dev.off()


